Structures and waves 
in a nonlinear heat-conducting medium 



Stefka Dimova, 

Faculty of Mathematics and Informatics, University of Sofia, Bulgaria 

dimova@f mi .uni-sof ia.bg 
Milena Dimova, Daniela Vasileva, 
Institute of Mathematics and Informatics, Bulgarian Acad. Sci. 
mkoleva, vasileva@math.bas .bg 

oo 

The final publication will appear in Springer Proceedings in Mathematics and 
Statistics Qhttp : //www . springer . com/ series/ 10533[ ) , Numerical Methods 
for PDEs: Theory, Algorithms and their Applications. 

<S . 

Abstract 

The paper is an overview of the main contributions of a Bulgarian team of researchers to the 
problem of finding the possible structures and waves in the open nonlinear heat conducting medium, 
described by a reaction-diffusion equation. Being posed and actively worked out by the Russian 
fT} | school of A. A. Samarskii and S.P. Kurdyumov since the seventies of the last century, this problem 

still contains open and challenging questions. 



O 

1 Introduction 



A very general form of the model of heat structures reads as follows: 

N 



ut = + Q( u )i t > °> x€R N , (1) 



i=i 



where the heat conductivity coefficients ki{u) > and the heat source Q(u) > are nonlinear functions 
of the temperature u(t, x) > 0. 

Models such as (JTJ) are studied by many researchers in various contexts. A part of this research is 
devoted to semilinear equations: fcj(ti) = 1, Q(u) = Xe u (Prank-Kamenetskii equation) and Q(u) = u^ 3 , 
f3 > 1. After the pioneer work of Fujita [28], these equations and some generalizations of theirs are 
studied intensively by many authors, including J. Bebernes, A. Bressan, H. Brezis, D. Eberly, A. 
Friedman, VA. Galaktionov, I.M. Gelfand, MA. Herrero, R. Kohn, LA. Lepin, SA. Posashkov, AA. 
Samarskii, J.L. Vazquez, J. J.L. Velazquez. The book [6] contains a part of these investigations and a 
large bibliography. 

The quasilinear equation is studied by D.G. Aronson, A. Friedman, HA. Levine, S. Kaplan, LA. 
Peletier, J.L. Vazquez and others. The contributions of the Russian school are significant. The unusual 
localization effect of the blow-up boundary regimes is discovered by numerical experiment in the work 
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of A. A. Samarskii and M.I. Sobol in 1963 [51]. The problem of localization for quasilinear equations 
with a source is posed by S.P. Kurdyumov [42] in 1974. The works of I.M. Gelfand, A.S. Kalashnikov, 
the scientists of the school of A. A. Samarskii and S.P. Kurdyumov are devoted to the challenging 
physical and mathematical problems, related with this model and its generalizations. Among them 
are: localization in space of the process of burning, different types of blow-up, arising of structures - 
traveling and standing waves, complex structures with varying degrees of symmetry. The combination 
of the computational experiment with the progress in the qualitative and analytical methods of the 
theory of ordinary and partial differential equations, the Lie and the Lie-Backlund group theory, has 
been crucial for the success of these investigations. The book [50j contains many of these results, 
achieved to 1986, in the review [M] there are citations of later works. 

A special part of these investigations is devoted to finding and studying different kinds of self- 
similar and invariant solutions of equation ([1]) with power nonlinearities : 

ki(u)=u a \ Q{u) = u p . (2) 

This choice is suggested by the following reasoning. 

First, such temperature dependencies are usual for many real processes [5], |54| . [57] . For example, 
when o"j = a = 2.5, j3 < 5.2, equation ([1]) describes thermo-nuclear combustion in plasma in the case of 
electron heat-conductivity; the parameters a = 0, 2</3<3 correspond to the models of autocatalytic 
processes with diffusion in the chemical reactors; <7~6.5 corresponds to the radiation heat-conductivity 
of the high-temperature plasma in the stars, and so on. 

Second, it is shown in [25], that in the class of power functions the symmetry of equation (fTJ 
is maximal in some sense - the equation admits a rich variety of invariant solutions. In general, 
almost all of the dissipative structures known so far are invariant or partially invariant solutions of 
nonlinear equations. The investigations of the dissipative structures provide reasons to believe that 
the invariant solutions describe the attr actors of the dissipative structures' evolution and thus they 
characterize important internal properties of the nonlinear dissipative medium. 

Third, this rich set of invariant solutions of equation ([I]) with power nonlinearities is necessary 
for the successful application of the methods for investigating the same equation in the case of more 
general dependencies ki(u), Q(u). By using the methods of operator comparison [29] and stationary 
states [32] it is possible to analyze the properties of the solutions (such as localization, blow-up, 
asymptotic behavior) of whole classes general nonlinear equations. The method of approximate self- 
similar solutions |31| . developed in the works of A. A. Samarskii and V.A. Galaktionov, makes it 
possible to put in accordance with such general equations some other, basic equations. The latter 
could have invariant solutions even if the original equations do not have such. Moreover, the original 
equations may significantly differ from the basic equations, and nevertheless their solutions tend to 
the invariant solutions of the basic equations at the asymptotic stage. 

Finally, in the case of power coefficients the dissipation and the source are coordinated so that 
complex structures arise, moreover, a spectrum of structures, burning consistently, occurs. 

Below we report about the main contributions of the Bulgarian research team: S.N. Dimova, 
IM.S. Kaschievl . M.G. Koleva, D.P. Vasileva, T.P. Chernogorova, to the problem of finding the possible 
evolution patterns in the heat-conducting medium, described by the reaction-diffusion equation ([1]), 
([2]). The outline of the paper is as follows. The main notions, needed further, are introduced in 
Section 2 on the simplest and the most studied radially symmetric case. The specific peculiarities 
of the numerical methods, developed and applied to solve the described problems, are systematized 
in Section 3. A brief report on the main achievements of the team is made in Section 4. Section 5 
contains some open problems. 
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2 The radially symmetric case, the main notions 

Let us introduce the main notions to be used further on the Cauchy problem for equation (pQ) with 
initial data: 

u(0,x) = uq(x) > 0, x G K , suptio(x) < oo. 

This problem could have global or blow-up solutions. The global in time solution is defined and 
bounded in l w for every t. The unbounded (blow-up) solution is defined in on a finite interval 
[0, To), moreover 

lim. t _^ T - sup u(t,x) = +oo. 



x 



The time To is called blow-up time. 

The unbounded solution of the Cauchy problem with finite support initial data uq(x) is called 
localized (in a strong sense), if the set 

U L = {ieR iV : u(T Q ,x) :=Iim"^ T - u(t,x) > 0} 

is bounded in R . The set Ql is called localization region. The solution localized in a strong sense 
grows infinitely for t — > T _ in a finite region 

u)L = {x G M. N : u(Tq,x) = oo} 



in general different from £1^. 
If for ki(u) = the condition 

r °° du 



l Q(u) 

holds, then the solution of the Cauchy problem is unbounded [50J. The heating of the medium happens 
in a blow-up regime, moreover the blow-up time of every point of the medium is different, depending 
on its initial temperature. 

If for Q(u) = the condition 

^-^du<+oo, i = l,2,...,N, (4) 
o u 

holds, then a finite speed of heat propagation takes place for a finite support initial perturbation in 
an absolutely cold medium [50] . 

In the case ([2]) of power nonlinearities it is sufficient to have o~i > 0, j3 > 1 for the conditions (|3|) 
and (|4|) to be satisfied. Then fej(0) = and equation ([T]) degenerates. In general it has a generalized 
solution, which could have discontinuous derivatives on the surface of degeneration {u = 0}. 

2.1 The basic blow-up regimes will be explained on the radially symmetric version of the Cauchy 
problem for equation (pQ): 

ut = ^ zl {x N - 1 u' T u x ) x + uP, xeR\, t>0, <r>0,p>l, (5) 

u x (t, 0) = 0, u(0, x) = u (x) > 0, < x < I, u (x) = 0, x > I. (6) 

If Uq(x) satisfies the additional conditions uq(x) G C(R^j_), (uqW o )(0) = 0, there exists unique local (in 
time) generalized solution u = u(t,x) of problem ([5j)-([6j) , which is a nonnegative continuous function 
in R\ x (0, T), where T G (0, oo] is the finite or infinite time of existence of the solution (see the 



3 



bibliography in the review [39]). Moreover u(t,x) is a classical solution in a vicinity of every point 
(t,x), where u{t,x) is strictly positive. It could not have the necessary smoothness at the points of 
degeneracy, but the heat flux — x N ~ 1 u a u x must be continuous. It means that u a u x = everywhere 
u = 0. 

Equation ([5]) admits a self-similar solution (s.-s.s.) [50j: 

u s (t, x) = <p(t)9 a ® = { 1 ~jr) ^ 0s(O, (7) 

z =x/m=x/ (i-±y-\ m= p-°-\ (8) 

The s.-s.s. corresponds to initial data u s (0,x) = 9 s (x). The function ip(t) determines the amplitude 
of the solution. The self-similar function (s.-s.f.) 9 S (£) > determines the space-time structure of 
the s.-s.s. ([7]). This function satisfies the degenerate ordinary differential equation in K+: 

m) s -^«"- W + + - f - o (») 

and the boundary conditions: 

#(0) = 0, ^(oo) = 0, W(£o) = 0, if S (£ O ) = 0. (10) 

-i 

Equation has two constant solutions: S (£) = 9h = (To(/3 — l))' 3-1 and a (£) = 0. These two 
solutions play an important role in the analysis of the different solutions of equation ©. For blow-up 
regimes we assume Tq > 0. Without loss of generality we set 

T = 1/09-1), then 9 H = 1. (11) 

The analysis of the solutions of problem © , ([TO]) , carried out in the works [39] , [52] , [26] , [1] , (see 
also [50J, Chapter IV), gives the following results: 

• For arbitrary 1 < /3 < a + 1 there exist a finite support solution 9 S (£) > 0. 

• For /3<<7 + l, iV > 1 and j3 = a + 1, N > 1 the problem has no nonmonotone solutions. The 
uniqueness is proved only for /3 < a + 1, iV = 1. 

The graphs of the s.-s.f. # s (£) f° r /3 = c + 1 = 3, iV = 1, 2, 3 are shown in Fig. 1, the graphs of 
the s.-s.f. 8 (£) for = 2.4 < a + 1 = 3, N = 1, 2, 3 - in Fig. 2. 

• For P > a + 1, N > 1 the problem has no finite support solutions. 

• If cr + 1 < /3 < /3 S = (<r + l)(iV + 2)/(iV — 2) + , f/3 s - i/ie critical Sobolev exponent), the problem 
has at least one solution 9 S (£) > in M 1 ^, strictly monotone decreasing in £ and having the 
asymptotics 

9 s (0 = C s C 2/if3 - a - 1) [l + u(0], w(£)->0, e^oo, (12) 
C s = C s (o~, f3, N) is a constant. Later on in [33] the interval in f3 has been extended. 

• For N = 1, /? > c + 1 the problem has at least 

* = -[-a]-l, a = / -1 - >1 (13) 
p — a — 1 

different solutions ([JJ, [26], [50J). Let us introduce the notations s ,i(€), i = 1,2,..., K for 
them. On the basis of linear analysis and some numerical results in the works [431 . |44j it has 
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2, P=2-4 




Figure 1: (3 — a + 1, 5- regime 



Figure 2: /3 < a + 1, i^S'-regime 



been supposed that the number of different solutions s ,i(£,) for /3 > a + 1 and N > 1 is -fT + 1. 
For N = 1 this result was refined [45 J by using bifurcation analysis: the number of solutions is 
K = [a], if a is not an integer, and K = a — 1, if o is an integer. For N = 2,3 the bifurcation 
analysis gives the same estimate for the number of different solutions, but for f3 ~ cr + 1, f3 > cr + 1 
it is violated (see 4.1). 



1.2 -I 




Figure 3: j5 > a + 1, LS'-regime 



Figure 4: /? > a + 1, LS'-regime 



The graphs of the four self-similar functions, existing for a = 2,0 = 3.6 > cr + 1 (K = 4) are shown 
in Fig. 3 (N = 1) and Fig. 4 (N = 3). 

These results determine the basic regimes of burning of the medium, described by the s.-s.s. 
0, ([HI). The following notions are useful for their characterization: 

- semi-width x s = x s (t), determined by the equation u(t,x s ) = u(t,0)/2 for solutions, monotone 
in x and having a single maximum at the point x = 0; 

- front-point Xf: u(t,xj) = 0, u a u x (t, Xf) = 0. 

2.1.1 ii^-evolution, total blow-up, 1 < /3 < a + 1 

The heat diffusion is more intensive than the heat source. The semi-width and the front tend to 



5 



infinity; a heat wave, which covers the whole space for time To, is formed. The process is not localized: 
mes = mes co L = oo, x s — > oo, xj — > oo, t — > T ~. 

2.1.2 5-evolution, regional blow-up, /3 = a + 1 

The heat diffusion and the source are correlated in such a way, that leads to localization of the 
process in a region Qi = ujl = {\x\ < L s /2} of diameter L s , called a fundamental length of the S- 
regime. The semi- width is constant, inside Ql the medium is heated to infinite temperature for time 
To : x s = const, Xf = L s /2. In the case N = 1 the solution O s (0 {Zmitrenko-Kurdyumov solution) is 
found [52J explicitly: 



2(a + l) 2 nC\ 1/a L a 
cos — , |£| < — 



a + 2 LsJ ~ 2 (14) 

|(| > ^, 



L s = diam 0,^ = (2iry/cr + 1)/<t, x s = L s arccos((2~ 2 )/tt). 

The solution (]14p is called elementary solution of the S-regime for iV = 1. In this case equation 
([9]) is autonomous and every function, consisting of k elementary solutions, k = 1,2,..., is a solution 
as well, i.e., equation ([9]) has a countable set of solutions. 

2.1.3 LS-evolution, single point blow-up, o~ + l<f3<f3f = a + l + jj 

Here /3/ is the critical Fujita exponent [46j . The intensity of the source is bigger, than the diffusion. 
The front of the s.-s.s. is at infinity (|12p . the semi- width decreases and the medium is heated to 
infinite temperature in a single point: 



mes ojl = 0, x s — > 0, t — > T ( 







According to the different s.-s.f. 9 St i(^), i = 1,2, ... , the medium burns as a simple structure (i = 1) 
and as complex structures (i > 1) with the same blow-up time. 

2.2 Stability of the self-similar solutions 

To show the important property of the s.-s.s. as attractors of wide classes of other solutions of the 
same equation, we will need of additional notions. 

In the case of arbitrary finite support initial data ^0(2;) © the so called self-similar representation 
[26J of the solution u(t,x) of problem (|5j), © is defined. It is determined at every time t according to 
the structure of the s.-s.s. ([7]), ([8]): 

6(U) = (1 - t/T )^ u - t/T )^) = ^-\t)u{t,^{t))- (15) 

The s.-s.s. u s (t, x) is called asymptotically stable [50], if there exists a sufficiently large class of solutions 
u(t,x) of problem ([5]), © for initial data uq(x) ^ 9 s (x), whose self-similar representations G(t,() tend 
in some norm to 6 S (£,) when t — > T ~~: 

||G(i,£)-0 s (£)||^O, t^T . (16) 

The definition of the self-similar representation (|15p contains the blow-up time To. For theoretical 
investigations this is natural, but for numerical investigations definition (|15p is unusable since for 
arbitrary initial data ^0(2;) 7^ 9 s (x) Tq is not known. Therefor another approach has been proposed 
and numerically implemented (for N = 1) in |49j . |26j . This approach gives a possibility to investigate 
the structural stability of the unbounded solutions in a special "self-similar" norm, consistent 
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for every t with the geometric form of the solution and not using explicitly the blow-up time To. A 
new self-similar representation, consistent with the structure of the s.-s.s. (J7j), ([8j) is introduced: 



eW) = .(«hMn/7( t ), *> = =^|f < 17 > 

If the limit (|16p takes place for 0(t,£), given in (|17p . then the self-similar solution u s (t,x) is called 
structurally stable. 

The notion of structural stability, i.e., the preservation in time of some characteristics of the 
structures, such as geometric form, rate of growth, localization in space, is tightly connected with the 
notion invariance of the solutions with respect to the transformations, involving the time [30]. This 
determines its advisability for investigating the asymptotic behavior of the blow-up solutions. 

In the case of complex structures another notion of stability is needed, namely metastability. 
The self-similar solution u s (t,x) is called metastable, if for every e > there exists a class of initial 
data u(Q, x) ~ s {x) and a time T, To — T <C To such that 

|| e(t,f)-0«(O ll< £ > for 0<t<T 

holds for the self-similar representations (|17p of the corresponding solutions. This means, that the 
metastable s.-s.s. preserves its complex space-time structure during the evolution up to time T, very 
close to the blow-up time To. After that time the complex structure could degenerate into one or 
several simple structures. 



3 Numerical methods 

To solve the reaction-diffusion problem ©, © and the corresponding self-similar problem @, (fTUj) . 
as well as their generalizations both for systems of such equations and for the 2D case, appropriate 
numerical methods and algorithms were developed. 

The difficulties, common for the nonstationary and for the self-similar problems were: the nonlin- 
earity, the dependence on a number of parameters (not less than 3); insufficient smoothness of the 
solutions on the degeneration surface, where the solutions vanish. In the case of radial symmetry for 
N > 1 and polar coordinates in the 2D case additional singularity at x = (£ = 0) occurs. 

The main challenge in solving the self-similar problems is the non-uniqueness of their solutions 
for some ranges of the parameters. The following problems arise: to find a "good" approximation to 
each of the solutions; to construct an iteration process, converging fast to the desired solution (corre- 
sponding to the initial approximation) and ensuring sufficient accuracy; to construct a computational 
process, which enables finding all different solutions for given parameters (a, j3, N) in one and the same 
way; to determine in advance where to translate the boundary conditions from infinity, for example 
(HDD, for the asymptotics (HD to be fulfilled. 

The main difficulty in solving the nonstationary problems is the blow-up of their solutions: blow- 
up in a single point, in a finite region and in the whole space. And two others, related with it: the 
moving front of the solution, where it is often not sufficiently smooth; the instability of the blow-up 
solutions. 

3.1 Initial approximations to the different s.-s.f. for a given set of parameters 

To overcome the difficulty with the initial approximations to the different s.-s.f., we have used the 
approach, proposed and used in the works [33], [Sj. Based on the hypothesis that in the region of 
their nonmonotonicity the s.-s.f. have small oscillations around the homogenious solution 8h, this 
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approach consists of "linearization" of the self-similar equation around Oh and followed by "sewing" 
the solutions of the resulting linear equation with the known asymptotics at infinity, e.g. (I12p . 

Our experiments showed [21], that when f3 — > a + 1 + the hypothesis about small oscillations 
of the s.-s.f. around Oh, is not fulfilled. The detailed analytical and numerical investigations [21J, 
[38J, [12J of the "linear approximations" in the radially symmetric case for iV = 1,2,3, showed that 
even in the case these approximations take negative values in vicinity of the origin, they still give the 
true number of crossings with 9h and thus, the character of nonmonotonicity of the s.-s. functions. 
Recommendations of how to use the linear approximations in these cases are made in [21] . 

Let us note, the "linear approximations" are expressed by different special functions: the confluent 
hypergeometric function iiq(a, b; z) and the Bessel function Jk(z) for different parameter ranges within 
the complex plane for a and b and different ranges of the variable z. To compute these special functions, 
various methods were used: Taylor series expansions, expansions in ascending series of Chebyshev 
polynomials, rational approximations, asymptotic series [22], [15] . 

3.2 Numerical method for the self-similar problems 

To solve the self-similar problem ([9]), (|10p and its generalization for systems of ODE and for the 2D- 
case, the Continuous Analog of the Newton's Method (CANM) was used ([H]-[2T], [38], [3D], [UJ). 
Proposed by Gavurin [35], this method was further developed in [56], [48] and used for solving many 
nonlinear problems. The idea behind it is to reduce the stationary problem L(6) = to the evolution 
one: 



by introducing a continuous parameter t, < t < oo, on which the unknown solution depends: 
8 = 0(£,t). By setting v = d9/dt and applying the Euler's method to the Cauchy problem (|18p . one 
comes to the iteration scheme: 



The linear equations (|19p (or the system of such equations in the case of a two-component medium) 
are solved by the Galerkin Finite Element Method (GFEM) at every iteration step. The combination 
of the CANM and the GFEM turned out to be very successful. The linear system of the FEM with 
nonsymmetric matrix is solved by using the LU decomposition. The iteration process (I20p converges 
very fast - usually less than 15-16 iterations are sufficient for stop-criterium ||L(0 n )|| < 10~ 7 . The 
numerical investigation of the accuracy of the method being implemented shows errors (i) of order 
0(/i 4 ) when using quadratic elements in the radially symmetric case, and (ii) of optimal order 0{h?) 
when using linear elements in the same case or bilinear ones in the 2D case. To achieve the same 
accuracy in vicinity of the origin in the radially symmetric case for N > 3 the nonsymmetric Galerkin 
method [27] was developed [38], [13] . 

The computing of the solutions of the linearized self-similar equation and their sewing with the 
known asymptotics is implemented in a software, so the process is fully automatized. The software 
enables the computing of the self-similar functions for all of the blow-up regimes, moreover in the case 
of LS-regime only the number k of the self-similar function O s ,k(Q must be given. 

3.3 Numerical method for the reaction-diffusion problems 

The Galerkin Finite Element Method (GFEM), based on the Kirchhoff transformation of the nonlinear 
heat-conductivity coefficient: 




L{d), 0(£,O) = O (O 



(18) 



L\e n )v n = -L(9 n ), 

0„+l = n + T n V n , < T n < 1, 72 = 0,1,..., 
On = On(0 = 0(C,t n ), V n = V n (C) = v($,t n ), 

#o(0 being the initial approximation. 



(19) 



(20) 




(21) 
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was used ([3], @], [22], [I2]-[l6], p2]-[20], [23], [53]) for solving the reaction-diffusion problems. This 
transformation is crucial for the further interpolation of the nonlinear coefficients on the basis of the 
finite element space and for optimizing of the computational process. 

Here below we point out the main steps of the method on the problem ([5]), (|6|) in a finite interval 
[0, X(t)] under the boundary condition u(X(t)) = 0. Because of the finite speed of heat propagation 
we choose X = X{t) so as to avoid the influence of the boundary condition on the solution. The 
discretization is made on the Galerkin form of the problem: 

Find a function u(t, x) G D, 

D = {u: x^-^u, x^ N -^/ 2 du^/ 2 /dx G L 2 , u(X{t)) = 0}, 
which for every fixed t satisfies the integral identity 

(u t ,v) = A(t;u,v), Vi> G -ff 1 (0, X(t)), 0< i < T , (22) 



and the initial condition 
Here 



(u,v) 



x(t) 



x N 1 u(x)v(x)dx, A(t;u,v) 



x(t) 



N _ x dG{u) dv a ' 
x — - — — h xvTv 



dx dx 



dx, 



H\0,X{t)) = {v : x^-^^v, x^-^v' G L 2 (0,X(t)), v(X(t)) = 0.} 



The lumped mass finite element method [55J with interpolation of the nonlinear coefficients G(u) f)21|) 
and q(u) = u 13 : 

n n 

G(u) ~ Gi = G(ui)ipi(x), q{u) ~ qi = q(ui)ipi(x) 

i=l i=l 

on the basis {<fi}, i = 1, . . . , n, of the finite element space is used for discretization of (|22[) . The result- 
ing system of ordinary differential equations with respect to the vector U (t) = (ui(t),U2(t), . . . , u n (t)) T 
of the nodal values of the solution u(t, x) at time t is: 



U = M~ 



-KG(U)) + q(U), U(0) = U . 



(23) 



Here the following denotations are used: G(U) = (G(u\), . . . ,G(u n )) T , q(U) = (q(ui), . . . ,q(u n )) T , 
M is the lumped mass matrix, K is the stiffness matrix. Let us mention, thanks to the Kirchhoff 
transformation and the interpolation of the nonlinear coefficients only the two vectors G(U) and q(U) 
contain the nonlinearity of the problem, while the matrix K does not depend on the unknown solution. 

To solve the system (f2"3"|) an explicit Runge-Kutta method [17] of second order of accuracy and an 
extended region of stability was used. A special algorithm for choosing the time-step r ensures the 
validity of the weak maximum principle and, in the case of smooth solutions, the achievement of a 
given accuracy e up to the end of the time interval. The stop criterion is r < f0~ 16 and then Tq is the 
approximate blow-up time, found in the computations. 

It is worth mentioning, that the nonlinearity has changed the prevailing opinion about the explicit 
methods. Indeed, there are at least two reasons an explicit method to be preferred over the implicit 
one for solving the system ([23]): 

- the condition for solvability of the nonlinear discrete system on the upper time level imposes the 
same restriction on the relation "time step - step in space", as does the condition for validity of the 
weak maximum principle for the explicit scheme (see [50J, chapter VII, §5); 

- the explicit method for solving large discrete systems has a significant advantage over the implicit 
one with respect to the computational complexity. 
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Let us also mention, that in the case of blow-up solutions the discrete system on the upper time 
level would connect solution values differing by 6-12 orders of magnitude, which causes additional 
difficulties to overcome. Finally the explicit methods allow easy parallelization. 

The special achievement of the proposed methods are the adaptive meshes in the LS-regime 
(refinement of the mesh) and in the H S-regime (stretching meshes with constant number of mesh- 
points), consistent with the self-similar low. Let us briefly describe this adaptation idea on the 
differential problem 

u t = Lu, u = u(t,x), i£R N , t>0, (24) 
which admits a self-similar solution of the kind 



u s (t, x) = ip(t)9s(£), Z = x/il>(t), u s (0, x) = e s {x). (25) 

Since the invariant solution u s (x,t) is an attractor of the solutions of equation (|24|) for large classes 
of initial data different from 9 s (x), it is important to incorporate the structure (|25p in the numerical 
method for solving equation (|24p . The relation (|25[) between £ and x gives the idea how to adapt the 
mesh in space. Let Ax^ be the initial step in space, Ax^ - the step in space at t = t k . Then Ax^ 
must be chosen so that A£( fe ) is bounded from below and from above 

Ax (0) /A < A^ k) < XAx^ 

for an appropriate A (usually A = 2). 

Further, by using the relation between ip(t) and (p(t), it is possible to incorporate the structure 
(|25p of the s.-s.s. in the adaptation procedure. In the case of equation © we have 

£ = xT(t) m , AZ = AxT(t) m , m = (f3 — <r — l)/2, T(t) = max ^ u ^ x ) . (26) 

max x uo(x) 



On the basis of the relations (126p the following strategy is accepted. 

In the case of a single point blow-up, m > 0, we choose the step Ax^ so that the step A£( fc ) be 
bounded from above: 

A^ fe ) = Ax {k) T(t) m < AAx (0) . (27) 

When condition (|27|) is violated, the following procedure is carried out: every element in the region, 
in which the solution is not established with a given accuracy 5 U (usually 5 U = 10 -7 ), is divided into 
two equal elements and the values of the solution in the new mesh points are found by interpolating 
from the old values; the elements, in which the solution is established with a given accuracy S u , are 
neglected. 

In the case of a total blow-up, m < 0, we choose the step Ax^ so that the step A^ k ^ be bounded 
from below: 

A£( fc) = Ax {k) T{t) m > AAs<°>. (28) 



When condition (J28J) is violated, the lengths of the elements are doubled, and so is the interval in x: 
X(tk+i) = 2X(tk), thus the number of mesh points remains constant. 

This adaptation procedure makes it possible to compute efficiently the single point blow-up as well 
as the total blow-up solutions up to amplitudes 10 6 — 10 12 , depending on the medium parameters. It 
ensures the authenticity of the results of investigation of the structural stability and the metastability of 
the self-similar solutions. Let us note, this approach does not require an auxiliary differential problem 
for the mesh to be solved, unlike the moving mesh methods [9j. The idea to use the invariant properties 
of the differential equations and their solutions [7j and to incorporate the structural properties (e.g. 
geometry, different kind of symmetries, the conservation laws) of the continuous problems in the 
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Figure 5: Graphs of the s.-s.f. S , 2 for N = 3, a = 2, = {3.6(1); 3.38(2); 3.2(3); 3.08(4); 3.03(5)} 
(left) and 6 S , 3 for N = 3, a = 2, /3 = {3.2(1); 3.1(2); 3.08(3); 3.06(4); 3.03(5)} (right). 



numerical method, lies at the basis of an important direction of the computational mathematics - 
geometric integration, to which many works and monographs are devoted - see [10], [21], [37] and the 
references therein. 

The numerous computational experiments carried out with the exact self-similar initial data fll4f> 
for N = 1, /3 = a + 1, as well as with the computed self-similar initial data, show a good blow-up time 
restoration (set into the self-similar problem) in the process of solving the reaction-diffusion problem. 
The preservation of the self-similarity and the restoration of the blow-up time demonstrate the high 
quality of the numerical methods for solving both the self-similar and the nonstationary problems. 

4 Results and achievements 

The developed numerical technique was used to analyze and solve a number of open problems. Below 
we present briefly some of them. 

4.1 The transition LS- to S-regime in the radially symmetric case 

The investigation of the limit case j3 — > cr + 1 + resolved the following paradox for N > 1: there exists 
one simple-structure s.-s. function in S'-regime (/3 = a + 1), whereas in LS-regime for /3 — > a + 1 + 
their number tends to infinity according to formula (|13|) . The detailed numerical experiment in [38J, 
|21j yielded the following results. First, it was shown, that the structure of the s.-s.f. for N > 1 and 
(3 ~ cr + 1, (3 > a + 1 is substantially different from the one for N = 1. Second, for N = 1 the transition 
(3 — )• <r + 1 + is "continuous" - the self-similar function s ,k(O f° r the LS-regime tends to a s.-s.f. of 
the /S'-regime, consisting of k elementary solutions. For N > 1 the transition behaves very differently. 

For a fixed and (3 — > a + 1 + the central minimum of the "even" s.-s.f. ®2j \ 3 = 1>2, ... , 
decreases and surprisingly becomes zero for some (5 = f3*-(a,N) (Fig. 5, left). For f3 < f3*(a,N) all 
of the s.-s.f. have zero region around the center of symmetry, and the radius of this region tends to 
infinity for (3 — > a + 1 + 0. All of the maxima of ^^-(C) tend to the maximum of the s.-s.f. of the 

S-regime for the corresponding a and for iV = 1. Thus the s.-s.f. &g*2j(£), "going to infinity" when 
(3 — > a + 1 + 0, tends to a s.-s.f. of the S'-regime for N = 1 and the same a, consisting of j elementary 
solutions. 

For fixed a there exists such a value f3j*(a,N), that for a + 1 < f3 < f3** the "odd" s.-s.f. 
^s2j+l(£)> 3 = 1) 2, . . . split into two parts: a central one, tending to the s.-s.f. of the S-regime for the 
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same N, and second one, coinciding with the s.-s.f. 6 S 2j(Qi "g° m g to infinity" when /3 — > a + 1 + 
(Fig. 5, right). 

According to the described "scenario", when f3 — > a + 1 + only the first s.-s.f. of the LS'-regime 
remains, and it tends to the unique s.-s.f. of the .S-regime. 

As a result of this investigation new-structure s.-s.f. were found - s.-s.f. with a left front. 
The existence of such s.-s. functions was confirmed by an asymptotic analysis, i.e., the asymptotics 
in the neighborhood of the left front-point was found analytically [12]. This new type of solutions 
initiated investigations of other authors [36] , |45j by other methods (the method of dynamical analogy, 
bifurcation analysis). Their investigations confirmed our results. 

4.2 The asymptotic behavior of the blow-up solutions of problem ([5]), ([6]) beyond the 
critical Fujita exponent 

For (3 > /3j = a+l+2/N the problem ([5]), ([6]) could have blow-up- or global solutions depending on the 
initial data. For the s.-s. blow-up solution ([7]), ([8]) it holds u s (t,r) L±(M N ). The qualitative theory 
of nonstationary averaging "amplitude-semi- width" predicts a self-similar behavior of the amplitude of 
the blow-up solutions and a possible non-self-similar behavior of the semi- width [50J . A question was 
posed there: what kind of invariant or approximate s.-s.s. describes the asymptotic stage (t —> Tq~) 
of the blow-up process? 

The detailed numerical experiment carried out in [23], [53] showed that the s.-s.s. ([7]), ([8j), corre- 
sponding to O s ,i(0i i s structurally stable: all of the numerical experiments with finite-support initial 
data ©, ensuring blow-up, yield solutions tending to the self-similar one on the asymptotic stage. 

4.3 Asymptotically self-similar blow-up beyond some other critical exponents 

The numerical investigation of the blow-up processes in the radially symmetric case for high space- 
dimensions N was carried out in [13] . [14] . The aim was to check some hypotheses [33] about the 
solutions of the s.-s. problem Q, (fTOj) and about the asymptotic stability of the corresponding s.-s. 
solutions for parameters beyond the following critical exponents: 
j3 s = (a + l)(N + 2)/(N - 2), N > 3 (Sobolev's exponent); 
I3 U = (a + 1)(1 + 4/(A - 4 - N > 11; 

fa = 1 + 3(<7 + 1) + (a 2 (N - 10) 2 + 2cr(5a + 1)(N - 10) + 9(<r + I) 2 ) 1 / 2 /(N — 10), N > 11. 

Self-similar functions, monotone in space, were constructed numerically for all of these cases, thus 
confirming the hypotheses of their existence (not proved for /3 > f3 u ). It was also shown, that the 
corresponding s.-s.s. are structurally stable, thus confirming another hypothesis of |33j . Due to the 
strong singularity at the origin, the nonsymmetric Galerkin method and the special refinement of the 
finite element mesh |14] were crucial for the success of these investigations. 

4.4 Two-component nonlinear medium 

The methods, developed for the radially symmetric problems ([5]), ([6]) and Q, (fTOj) were generalized 
in [TU], |40j . [53| for the case of two-component nonlinear medium, described by the system: 

x£R l + , N = 1,2,3, 

(29) 

Gi > 0, ^ > 1, 7i > 0,i = 1,2. 
This system admits blow-up s.-s.s. of the form 

u u = (i - t/T r^e ls (o, z = x/(i - t/nr, 

(30) 

u 2s = (1 - t/ToHMO, mt < 0, i = 1, 2, 



uu = -^zj{x N l u a l 1 u lx ) x + u^uf, 

1 / N—l cro \ i 7i 02 
U 2t = -J^ii x U 2 2 U 2x ) x + U{ U'z , 
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where 

mi = — , ai = 7j + 1 - ft, i = 1, 2, p = (ft - l)(ft - 1) - 7i72, 
P 

mio"! + 1 m 2 a 2 + 1 , , , . 

n= 5 = 2 ' + 1 - ft) = ^2(71 + 1 - Pi)- 

The s.-s.f. satisfy the system of nonlinear ODE: 

1 
1 



HOu, o 2s ) = ~?h : (Z N - 1 o£e , uy + nte> ls - mi e ls - efretg = o, 
WOu, o 2s ) = -J^ie-^teLY + - m 2 e 2s - e%o% = o 



(31) 



and the boundary conditions 

hm = 0, lim 9 is = 0, z = 1,2. (32) 

Superconvergence of the FEM (of order (3(/i 4 )) for solving the s.-s. problem (I3ip -( [32"1) by means of 
quadratic elements and optimal-order convergence (0(h 2 )) by means of linear elements were achieved. 
The structural stability of the s.-s.s. (i30j) for parameters crj, ft, 7$, i = 1,2, corresponding to the 
LS'-regime (n > 0) was analyzed in |40j, [53], [19]. It was shown that only the s.-s.s. of systems (|29p 
with strong feedback (p < 0), corresponding to the s.-s.f. with two simple-structure components, were 
structurally stable. All the other s.-s.s. were metastable - self-similarity was preserved to times not 
less than 99.3%T . 

The proposed computational technique can be applied to investigating the self-organization pro- 
cesses in wide classes of nonlinear dissipative media described by nonlinear reaction-diffusion systems. 

4.5 Directed heat diffusion in a nonlinear anisotropic medium 

Historically the first Bulgarian contribution to the topic under consideration was the numerical real- 
ization of the self-similar solutions, describing directed heat diffusion and burning of a two-dimensional 
nonlinear anisotropic medium. It was shown in [25] that the model of heat structures in the anisotropic 

case 

u t = {u ai u Xl ) Xl + {u a2 u X2 ) X2 +u p , x = {xx,x 2 ) E M 2 , a x > 0, a 2 > 0, > 1 (33) 
admits invariant solutions of the kind 



u s (t,xi,x 2 ) = (1- — ) e s (0, £ = (6,6) e 



d2 



Xi/[i- — \ , Tm = - ,1 = 1,2. 



The self-similar function # s (£i>6) satisfies the nonlinear elliptic problem 

9 fna, 9d s\ , /3-aj-l ,de, 



m) = ±(-±( c %) + = 0, (34) 



i=i ^ 

d9 

— s - = 0, i = 1,2; e s {i) 0, |e| -> 00. (35) 

The Cauchy problem for equation (|33p was investigated in the works [3], [I] for different parameters 
o"i, o"2 and ft Depending on the parameters, different mixed regimes: S — HS, HS — LS, S — LS of 
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Figure 6: S'-regime: a± = 2, cr 2 = 2, (3 = 3 




E1M2 E1M3 E1M4 




E1M5 E1M5 E1M7 

Figure 9: Self-similar functions EjMm, j = 1, m = 2, 3, 4, 5, 6, 7, a = 2, P = 3.25. 



heat transfer and burning were implemented numerically. The evolution in time of one and the same 
initial perturbation is shown for the cases of the 2D radially symmetric S-regime (Fig. 6), the mixed 
HS — S-regime (Fig. 7) and the mixed HS — L S-regime (Fig. 8). 

To solve the Cauchy problem for equation (133p a modification of the TERMO Package of Applied 
Programs [8], designed initially for solving isotropic problems with piecewise constant coefficients, was 
done. TERMO had been worked out by an IMFBAS team, after the idea of Raytcho Lazarov and 
under his guidance, a merit worth mentioning here. 

The s.-s. functions for the corresponding mixed regimes were found in [3] by self-similar processing 
of the solution of the Cauchy problem for equation (|33l) . Later, in [IT], they were found as solutions of 
the self-similar problem (|34p . (|35p . The self-similar functions of complex symmetry for the isotropic 
case in Cartesian coordinates (denoted in [44] as Ei/j) were found as a special case. 

Later on, in |41| . |40j . the numerical methods were modified for the isotropic 2D self-similar 
problem in polar coordinates to construct numerically another class of self-similar functions of complex 
symmetry (denoted in [33] as EjMm) in LS'-regime and to investigate their structural stability. 

Graphical representations of the evolution of the anisotropic invariant solutions, as well as the s.-s. 
functions EjMm for some different values of j,m,a,f3, are included in the Handbook [11] . We show 
some of the s.-s.f. EjMm in Fig. 9. 

4.6 Spiral waves in //S'-regime 

The numerical realization of the invariant solutions, describing "spiral" propagation of the nonhomo- 
geneities in two-dimensional isotropic medium appears to be one of the most interesting contributions 
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of ours. The mathematical model in polar coordinates reads: 



u t = ^(ru a u r ) r + -^(ifu^^ + u p , a > 0, > 1. (36) 
It admits s.-s.s. of the kind [I], [30J: 

u s (t,r,ip)= (l-^) Bs{^), (37) 

S-rlU-lY*, i- v + J*-M(l-'), m-^^pi. (38) 



V r /3-i v V 2 

The self-similar function S (£,</>) satisfies the nonlinear elliptic equation 



e V 3£ / £ 2 5<a \ ' d<i>) 2 ? d£ 



(39) 



0-1 

Here cq 7^ is the parameter of the family of solutions. From (|38|) it follows 



sip 0-a-l 



£e s<p = re Sip = const, s 



2c 



This means, that the trajectories of the nonhomogeneities in the medium (say local maxima) are 
logarithmic spirals for /3 7^ a + 1 or circles for = a + 1. The direction of movement for fixed cq, 
for example cq > 0, depends on the relation between a and : for > a + 1 - towards the center 
(twisting spirals), for < a + 1 - from the center (untwisting spirals). 

The problem for the numerical realization of the spiral s.-s.s. ()37|) . (|38p was posed in 1984, when 
the possibility for their existence has been established by the method of invariant group analysis in 
the PhD Thesis of S.R. Svirshchevskii. As it was stated in |2j. there were significant difficulties for 
finding such solutions. First, the linearization of the self-similar equation (|39p was not expected to 
give the desired result, because it is not possible to separate the variables in the linearized equation. 
Second, the asymptotics at infinity of the solutions of the self-similar equation were not known. 

The first successful step was the appropriate (complex) separation of variables in the linearized 
equation. Using the assumption for small oscillations of the s.-s.f. S (£,0) around 9 l H = 1, i. e., 
0s(£ 5 0) = 1 + a u{£,i </>)> a = const, I ay I <C 1 and the idea of linearization around it, the following 
linear equation for y(£,(j>) was found [22j: 

1 d ( dy\ 1 d 2 y — a — 1 By dy . . 

~m ) ~ eW + - co ^ + (1 - p)y = °- 

Seeking for particular solutions Yfc(£, 4>) = R^^e 1 ^ , k € N, it was found: for = a + 1 

where Jk(z) is the first kind Bessel function of order k, and for ^ a + 1 

Rk(0=Z k lF 1 (a,b;z), a = -t— 5L- + -, 6 = 1 + fc, ^ = ^ f, 

p — a — 1 I 4 
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k = 1 k = 2 k = 3 

Figure 10: One-armed spiral solution {k = 1), two-armed spiral solution (k = 2), three-armed 
spiral solution (k = 3), cq = 1, a = 3, /3 = 3.6. 



where iiq(a,&, z) is the confluent hypergeometric function. The detailed analytical and numerical 
investigation [22] of the functions yfc(£,<^>) = 3f?(Yfc(£, 6)) showed that their asymptotics at infinity are 
self-similar, as well as that the functions 

9 s , k {£,<f>) = l + uy k (£,(f>), |oy fc |«l (40) 

are very close to the sought-after solutions 9 S (^,6). Moreover, the amplitude of the linear approxi- 
mations 2/fc(^,</ ) ) tends to zero for £ — > oo in the case of -ffS-regime, and to infinity in the case of 
LS-regime. This gave the idea to seek for s.-s.s. of the H ^-regime, tending to the nontrivial constant 
solution 9 S = Oh-, i-e., to generalize the notion of s.-s. functions, and consequently, the notion of the 
structures and waves, which arise and preserve themselves in the absolutely cold medium. Although 
not exploited earlier, this change is reasonable and fully adequate to the real systems. 

All stated above enabled us to predict the asymptotics of the solutions of equation (139f) . to derive a 
boundary condition by using this asymptotics and to close the s.-s. problem by the following boundary 
and periodic conditions [15] . [15] . [20] : 



hm^fe^- 
99 Sy k 



8,k 



7 S,fc 



= 0, ct>e 

1 ^k 



2vr 



sm 



+ -Int + fj,), £ = I > 1, <f>e 
s 



2vr 



(41) 



06 l ^' Uj dch I k ) 



o < i < /, 



where in = m/((3 — 1), 7 and \x are constants, depending on <r, /3,co- The numerical solving of the 
problem (j39h . (I4ip for Co 7^ with initial approximations (140j) gives "the spiral" s.-s. functions of the 
-ffS'-regime; some of them and their evolution in time by solving equation (|36p were investigated in 
[18J, [20], [15] . The graphs of the s.-s.f. for k = 1 (one-armed spiral), k = 2 (two-armed spiral), k = 3 
(three-armed spiral) are shown in Fig. [TU1 The rest of the parameters are a = 3, f3 = 3.6, cq = 1. 

The evolution of three-armed s.-s.f. for parameters a = 2, /3 = 2.4, /c = 3, Co = 1 is shown in 
Fig.HU The exact blow-up time is (llip To = l/(/3— 1) = 0.(714285). Similarly to all complex-symmetry 
s.s.-s., the three-armed spiral one is metastable - at T — > Tq it degenerates into the simplest radially 
symmetric s.-s.s. for the same parameters a, (3. The mesh adaptation in r-direction when solving 
equation (|36p is realized again in consistency with the self-similar law, keeping the same number 
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Figure 11: Evolution of three-armed spiral wave: a = 2, (3 = 2.4, c = 1, k = 3, T = 0.(714285). 



of mesh points during the whole process of evolution. Having to solve two nonlinear 2D problems 
(elliptic and reaction-diffusion), going through a number of approximations, it is astonishing that the 
restoration (with accuracy 10 -6 ) of the exact blow-up time is practically perfect (Fig. 11, right most 
bottom). 

The existence of spiral structures in LS'-regime is an open question by now. The "ridges" of the 
linear approximations of the s.-s.f. in the LS'-regime tend also to the self-similar ones for £ — > oo 
[22J, but their amplitudes tend to infinity and it is not clear by what asymptotics to sew the linear 
approximations. 

4.7 Complex symmetry waves in H S-regime and S-regime 

In the process of solving the problem for the spiral s.-s.s. an idea arose to seek for complex non- 
monotone waves, tending to the nonzero homogeneous solution for £ — > oo in iLS-regime and Co = 0. 
Fig. [12] shows a complex-symmetry s.-s.f. and its evolution in time for the same parameters a, f3, as in 
Fig. HH and k = 2, cq = 0. Note the same perfect restoration of the blow-up time. The results about 
the spiral and the complex-symmetry s.-s.f. are included in the book |54j . 

Later the problem of finding nonmonotone s.-s.f. in the S-regime, tending to the nontrivial ho- 
mogeneous solution Oh-, was posed and successfully solved [16]. The self-similar equation was solved 
with boundary conditions 
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The s.-s.f. for /3 = a + 1 = 3 and its evolution in time are shown in Fig. [T3j 
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Figure 12: Evolution of a complex wave in ifS-regime: a = 2, (3 = 2.4, Co = 0, k = 2, T 
0.(714285). 



Let us mention, that the existence of continuum of solutions to the radially symmetric s.-s. prob- 
lem in HS- and /S-regimes, which tend to the nontrivial homogenious solution Oh for £ — > oo, was 
mentioned in |50| . but this result has remained without attention. As it turned out, it is these solutions 
that determine the spiral structures and the complex structures in HS- and S-regimes. 



5 Open problems 

In spite of the numerous achievements related to the problems considered here, many interesting 
questions are still open: how many complex-symmetry s.-s.f. of the kind Ei/j and EiMj, tending to 
the trivial constant solution 9 S = 0, exist; how can the self-similar problem for the spiral s.-s.f. in the 
LS-regime be closed, and therefore how to construct these spiral s.-s.f. numerically; how wide are the 
different classes of spiral s.-s.f. for (3 < a + 1 and the different classes of complex-symmetry s.-s.f. in 
S- and if S'-regimes. 

Although the method of invariant group analysis shows that the differential problem admits some 
kind of self-similar solutions and their numerical realization is a constructive "proof of their existence, 
theoretical proofs are still missing in most of the cases. 

The numerical investigation of the accuracy of the approximate solutions on embedded grids shows 
optimal-order- and even super convergence results, but it would be interesting to find theoretical esti- 
mates as well. 

All these questions pose challenging problems both from theoretical and computational points of 
view. 
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Figure 13: Evolution of a complex wave in 5-regime: a = 2, j3 = 3, cq = 0, k = 2, Tq = 0.5. 
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